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Abstract: We update predictions for lepton fluxes from the hadropro duct ion of 
charm quarks in the scattering of primary cosmic rays with the Earth’s atmosphere. 
The calculation of charm-pair hadroproduction applies the latest results from pertur¬ 
bative QCD through next-to-next-to-leading order and modern parton distributions, 
together with estimates on various sources of uncertainties. Our predictions for the 
lepton fluxes turn ont to be compatible, within the uncertainty band, with recent 
resnlts in the literature. However, by taking into acconnt contributions neglected 
in previous works, our total uncertainties are mnch larger. The predictions are cru¬ 
cial for the interpretation of resnlts from nentrino experiments like IceCnbe, when 
disentangling signals of nentrinos of astrophysical origin from the atmospheric back¬ 
ground. 
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1 Introduction 

Atmospheric lepton fluxes are important backgrounds in the search of neutrinos of 
astrophysical origin [1]. In particular recent claims from the IceCube experiment, 
which has detected a statistically significant sample of leptonic events at very high 
energies [2, 3], whose interpretation is still under debate [4, 5], require an estimate 
of the background as accurately as possible. One of the most uncertain components 
of this background is the prompt contribution due to the hadroproduction of charm 
quarks in the hard scattering of primary cosmic rays with the Earth’s atmosphere, the 
so-called atmospheric charm. In this paper we will concentrate on the contribution 
to the lepton fluxes that can be ascribed to atmospheric charm. 

After initial studies on the basis of phenomenological models (see e.g. Ref. [6, 7] 
and references therein), previous predictions for lepton fluxes from atmospheric 
charm have been obtained within the framework of perturbative Quantum Chro¬ 
modynamics (QCD) for proton-proton collisions according to the standard QCD 
collinear factorization formalism, with hard-scattering evaluated at leading order 
(LO) in Ref. [8] and including radiative corrections at next-to-leading order (NLO) 
in Ref. [9], respectively. As an alternative description motivated by the high colli¬ 
sion energies of the underlying hard scattering. Ref. [10] has used the so-called color 
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dipole approach as an effective model for the production of colored particles at high 
energies, in order to compute the production rates for atmospheric charm. 

All these predictions, however, are subject to very large theoretical uncertainties. 
While the results of Ref. [10] are very sensitive to the parameters of the phenomeno¬ 
logical model for the color dipole, which are poorly constrained by experimental 
data, also the standard perturbative QCD predictions for the hadroproduction of 
charm quarks acquire big uncertainties, of the order of several ten percents, in the 
kinematical regions of interest for astrophysical applications. In the latter case, these 
uncertainties are due to estimates of the missing radiative corrections at higher or¬ 
ders, the knowledge of parton distribution functions (PDFs), especially the gluon 
PDF at small fractions x of the momenta of the colliding protons, as well as the 
precision on the charm quark mass. 

Since the start of the Large Hadron Collider (LHC) and thanks to both theo¬ 
retical and experimental progress, our understanding of charm-pair hadroproduction 
at high energies has significantly improved. It is, therefore, the aim of the present 
paper to provide new predictions for atmospheric charm and its contribution to lep¬ 
ton fluxes, on the basis of standard perturbative QCD, taking into account the most 
recent developments in this field. 

For inclusive charm-pair hadroproduction we use QCD predictions up to next- 
to-next-to-leading order (NNLO) in order to establish the apparent convergence of 
the perturbative expansion and the stability under variation of the renormalization 
and factorization scales, together with recent determinations of PDFs compatible 
with constraints from LHC measurements. We also investigate the dependence on 
the renormalization scheme and the value for the charm quark mass and discuss dif¬ 
ferences between the running mass and the pole mass schemes. A consistency check 
is performed by a comparison of the theory predictions to available LHC data from 
ALICE [11], ATLAS [12] and the LHCb [13] experiments obtained in the runs at 
y/S = 7 and 8 TeV center-of-mass energy, because the data are within the kinematic 
region of interest for atmospheric charm. The differential distributions for charmed 
hadron production which are necessary in order to compute the lepton fluxes are 
obtained in perturbative QCD with a consistent matching between NLO QCD cor¬ 
rections and parton showers, as the respective predictions at NNLO are currently not 
available. The interface to the PYTHIA event generator [14] accounts for the full effect 
of parton showers and the hadronization. Our study features a detailed discussion 
of the different sources of theoretical uncertainties affecting predictions, which are 
propagated to the computation of the lepton fluxes. In this way, a total uncertainty 
band for the final predictions of the prompt lepton fluxes is established and compared 
to previous results in the literature. The effects on the fluxes due to modifications 
in primary cosmic ray spectra are also shown by making use of the latest spectra 
available in an analytic form. 

Recently, the authors of Ref. [9] have proposed an update of that work in Ref. [15] . 
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Our computation is independent and differs from Ref. [15] because of the up-to-date 
perturbative QCD results and methods used in the computation of charm and D- 
hadron production cross-sections and because of the choice and the variation of the 
input parameters. In summary, this leads to a more comprehensive estimate of the 
related uncertainties for the prompt lepton fluxes. 

The paper is organized as follows: In Sec. 2 we present the method, the input, and 
the tools we have used for the calculation, together with examples of results from 
its intermediate steps. Sec. 3 contains our predictions for the lepton fluxes along 
with a discussion of the related uncertainties. In Sec. 4 we sketch the astrophysical 
implications of these predictions, after comparing them to those used so far by the 
astrophysical community, and discuss their potential implications for the IceCube 
experiment. Finally, in Sec. 5 we draw our conclusions, with reference to future 
theory progress and measurements which could help to decrease the uncertainties on 
the predictions presented. 


2 Method: cascade equations and their solution 


The particle evolution through an air column of depth X in the Earth’s atmosphere 
can be obtained by solving a set of coupled differential equations, so-called cascade 
equations. Following Ref. [9, 16] one has 


d(j)j 

dX 


\ ^ T ^ ^ Sprod{k } j) T ^ ^ Sdecay{k ^ j) T ^reg{j t j) . 

X,int X,dec 

( 2 . 1 ) 


A dependence on the energy Ej is understood in all terms of eq. (2.1), j labels a 
particle species, Xj^int and Xj^dec its interaction and decay lengths, respectively, while 
Sprod and Sdecay denote the generation functions for production and decay: 
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rsj 






(2.3) 


Here, (j)k{Ek,X) is the flux of particle Zc, ak is the total inelastic cross-section for the 
interaction of particle k in the atmosphere, dak^j/dEj is the energy distribution of 
particle j produced by fc, Fj is the total decay width of particle j and dVj^i/dEi is the 
energy distribution of particle I produced by the decay of j. Regeneration functions 
in eq. (2.1), i.e., SregU —t j), can be viewed as a particular case of Sprodik —>■ j) when 
k = j. According to the nature of particle j (nucleon, heavy-hadron, neutrino), some 
of the terms in eq. (2.1) may be absent. 
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The right hand sides of eqs. (2.2) and (2.3), dehning the so-called Z-moments 
for the prodnction and decay of particle j, respectively, are obtained after noticing 
that the X dependence of fluxes approximately factorizes from their E dependence. 
In this approximation, analytic solutions exist for eq. (2.1) in the limit where the 
energy of intermediate particles leading to hnal leptons is either very small or very 
large with respect to their critical energy, the latter being proportional to the par¬ 
ticle mass m and to the inverse of its proper life-time tq. In the vertical direction, 
Ecrit = m ho/{cTo), where ho is the vertical depth of the atmosphere, for which 
an isothermal model is assumed with the density of the atmosphere evolving with 
depth as p{h) = po exp(—h/ho). 

In fact, the competition between hadron interaction and hadron decay is crucial 
in determining the hnal lepton huxes and Ecru represents an approximate energy 
above which the hadron decay probabilities are suppressed with respect to their in¬ 
teraction probabilities. In particular, one can distinguish between the conventional 
neutrino hux and the prompt neutrino hux, according to the nature of the interme¬ 
diate hadrons. The conventional hux is produced by the decays of charged kaons 
and pions which dominate over their interaction rates at relatively low energies, as 
the critical energies for these particles are smaller than 1 TeV. On the other hand, 
for larger energies, the probability of secondary interactions overcomes the probabi¬ 
lity that these mesons decay, thereby progressively suppressing the hux of neutrinos 
from them. At energies above 10^ — 10® GeV, neutrinos are thus mainly produced 
by the decay of other particles. In the framework of the Standard Model, these are, 
in particular, charmed and bottomed heavy-hadrons, which are characterized by a 
larger critical energy {Ecru > 10^ GeV) than pions and kaons ^ . These immediately 
decaying particles (r ~ 10“^^ s) give rise to the so called prompt hux, that is the 
object of study of this paper. 

In case of hadrons decaying into leptons, the hux of leptons coming from low 
energy hadrons, i.e., from hadrons with E Ecrit^ can be approximated by 

(2.4) 

^pp 

whereas the hux of leptons from hadrons with E S> Ecrit is approximated by 

j^high _ ryhigh Eph Ecrit,h hl(A/i/Ap) ^^ 

— ^hl t 7 1 G-OJ 

^ ^pp -^h 1 — -T— 

with Aj{Ej) dehned as Aj{Ej) = \j{Ej)/{l — Zjj{Ej)). In the approximated solutions 
to the cascade equations outlined above, an energy dependence is understood in all 

^ More precisely, the critical energies in vertical direction for the charmed hadrons considered 
in this work amount to: = 9.71 -10^ GeV, = 3.84 -10^ GeV, = 8.40 -10^ GeV, 

E<^^t ^ 24.4 -10^ GeV. 
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fluxes (f), all Z-moments and all attenuation lengths A. Note that the low energy 
lepton flux is isotropic, whereas the high energy lepton flux is characterized by an 
angular dependence f{9) ~ 1/ cos(6*) for 9 < 60°, where 9 is the angle with respect 
to the zenith, and by a more complex angular dependence close to the horizon. 

The solution in the intermediate energy range E ~ Ecru is obtained by the 
geometric approximation 





( 2 . 6 ) 


whose quality and validity depend on the particular shapes of and as a 
function of the lepton energy, see, e.g.. Fig. 10 in Sec. 3 below, for an example of 
this interpolation. In this way one gets the contribution (j)h^i to the lepton flux from 
each hadron species h. Summation over all hadron species hnally provides the total 
lepton flux (pi for each lepton species I, that is (pi = Alternatively, the 

system of differential equations in eq. (2.1) can also be solved numerically. 


2.1 Input: cosmic ray primary spectrum 

The primary cosmic ray (CR) spectrum is an important input of our calculation as it 
enters the solution of the cascade equations (2.4) and (2.5) both implicitly through 
Z-moments and explicitly. 

CR spectra in the upper layer of our atmosphere are very well constrained 
at low energies by many direct measurements. Balloon-borne and space experi¬ 
ments, like AMS and CREAM, are able to discriminate with high precision the 
individual elements included in the cosmic ray composition up to energies around 
Eiab = 10® GeV [17]. On the other hand, the high energy tail of cosmic rays is 
subject to signihcant uncertainties, in particular related to the different possible 
CR compositions (protons or heavier ions, up to iron) and the CR origin (galactic or 
extra-galactic). The high energy region is investigated by ground-based experiments, 
like KASCADE, KASCADE-Grande and the Pierre Auger Observatory, which glo¬ 
bally cover the energy range between Eiab = 10® GeV and up to several 10^^ GeV. 
In this context it is important to note that the lepton fluxes at a given energy are 
affected by GR spectra at energies even larger by a factor of order 0(100 — 1000), 
due to the integration over primary energies in the expressions of the generation 
functions, eqs. (2.2) and (2.3). Therefore, in order to parametrize the uncertainty on 
our knowledge of cosmic ray spectra at high energies we consider the following pos¬ 
sibilities, i.e., we evaluate lepton fluxes separately for each of the following primary 
cosmic ray spectra, which are available in literature 

^Given the fact, that we parameterize the p-A cross-sections in terms of p-p ones, cf. Sec. 2.3, 
we consider the all-nucleon version of each spectrum. 
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1) Power-law spectrum, composed by two parts: 

= 1.7 E-2-^cm-2 s-^ sr-^ GeV^ for E < 5 ■ 10® GeV, 

174 s"^ sr"^ GeV^ for E > 5 ■ 10® GeV. (2.7) 

This is a reference spectrum used in earlier works on prompt lepton fluxes, and, 
although recent measurements have shown that it basically overestimates nu¬ 
cleon fluxes at the highest energies, we consider it for reference and comparison 
with older works [9, 10, 15]. 

2) Gaisser 2012 (variant 1 and 2) [18]: 

The hrst variant of the Gaisser spectrum, htting available experimental data 
of different origin to an analytic expression with a number of parameters, is 
based on the hypothesis that three populations, one including GR particles 
accelerated by SuperNova remnants in our galaxy, a second one still of galactic 
origin but with an higher energy, and a third one of particles accelerated at 
extra-galactic sources, contribute to the measured GR spectrum. The three 
populations are characterized by different rigidities they all include protons 
and nuclear groups (He, GNO, Mg-Si, Fe) with different spectral indices. The 
second variant of the Gaisser spectrum provides a special treatment of the 
third population, which is assumed to be composed of protons only, with large 
rigidity. 

3) Gaisser 2014 (variant 1 and 2) [19, 20]: 

This uses the same functional form as in Gaisser 2012, but with updated pa¬ 
rameters for an alternative £t of experimental data. In particular, the hrst 
variant of the spectrum involves three populations, two of galactic and one of 
extra-galactic origin, involving the p. He, G, O, Fe nuclear groups, with diffe¬ 
rent rigidities with respect to the Gaisser 2012 case. The second variant differs 
from the hrst one because it includes an additional component from heavier 
nuclei, plus a fourth population, characterized by extra-galactic protons only, 
with large rigidity. This ahects the ultra-high-energy part of the spectrum 
and improves the agreement with Auger data on cosmic ray composition at 
high-energy [21]. 

The all nucleon spectra corresponding to these cases are shown in Fig. 1. The ehect of 
the diherent options on the shape of lepton huxes is extensively discussed in Sec. 3.1. 
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Cosmic Ray primary aii-nucieon fiux 



Figure 1. The all-nucleon primary cosmic ray spectra used as input in this work. See 
text for more detail. 


2.2 Input: p-Air total inelastic cross-section 

The total inelastic proton-Air cross-section as a function of the laboratory energy 
is an input in the denominator of the integrand of the generation function Sprod 
in eq. (2.2), and as a consequence, for the Z-moments for heavy-hadron hadropro- 
duction. Several measurements exist for this quantity, performed by different ex¬ 
periments (for a collection of results see Ref. [22] and references therein), together 
with theoretical predictions, on the basis of phenomenological models. In this paper, 
for compatibility with previous works, we consider both the analytical formula [23], 
already used in old estimates of prompt neutrino fluxes (see, e.g., [9]), 

= 290 - 8.71n(Z/GeV) + l.U\n^{E/GeY) mb, (2.8) 

and predictions from the SIBYLL2.1 [24] and the QGSJetO. Ic [25] models for hadronic 
interactions included in the CORSIKA package [26]. Cross-sections corresponding to 

^The rigidity of each population multiplied by the atomic number of each nuclear group, de¬ 
termines the characteristic energy where the corresponding all-particle CR spectrum exponentially 
cuts off. The larger the rigidity is, more extended is the spectrum at high energy. 
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E|ab (GeV) 

Figure 2. The p-Aii total inelastic cross-section according to different models 
(QGSJetO. Ic, SIBYLL2.1, analytical) as compared to measurements from Auger and older 
experiments. The experimental data are taken from Ref. [22] and references therein. 


those possible options are shown in Fig. 2 together with presently available experi¬ 
mental data. We also point ont that predictions from other CORSIKA models, like 
EPOS 1.99 [27] lie within the band that one can draw from these two choices, as 
shown in Fig. 2 of Ref. [22], so that we consider them as npper and lower limits. 
The recent measnrement from the Anger collaboration at a/s = 57 TeV, reported in 
Ref. [22] tnrns ont to be in agreement, within the error bands, with the predictions 
from QGSJetO. Ic. 

Discnssions on the effects of the different assnmptions for the p-Air cross-section 
on onr hnal resnlts of lepton flnxes are reported in Sec. 3.2. 

2.3 Charm hadroproduction cross-section 

Heavy-qnark hadroprodnction has been extensively stndied in pertnrbative QCD. 
The QCD corrections at NLO have hrst been obtained in Refs. [28-30] and are avai¬ 
lable in pnblic tools, like hvqmnr [31], MCFM [32] or HELAC-NLO [33] for the antomatic 
compntation of fnlly differential observables. For the inclnsive cross-section, the 
QCD corrections are complete to NNLO [34-37] and, thns far, have been applied to 
top-qnark pair prodnction. All these theory predictions have adopted the on-shell 
renormalization scheme for the heavy-qnark mass. The conversion to the MS scheme 
for the heavy-qnark mass has been discnssed in Refs. [38-40]. 

Beyond the pertnrbative expansion at hxed order, the resnmmation of large lo¬ 
garithms featnres important improvements, cf., e.g., the review in Ref. [41] for charm 



















and bottom production at the LHC. In dynamical regimes where the transverse mo¬ 
mentum pt of the heavy quark is much larger than its mass m the semi-analytical 
resummmation of logarithms in pT/ra has been performed in Ref. [42] in the so- 
called FONLL approach. On the other hand, when the NLO corrections are consis¬ 
tently matched with parton showers (PS) as in the POWHEG [43, 44] or MC@NLO [45] 
approaches using the frameworks of POWHEG-BOX [46] or (a)MC@NL0 [47], respec¬ 
tively, the leading logarithms are effectively resummed through the Monte Carlo 
(NLO -P PS) event generators. 

In summary, there exists a robust theoretical framework with a set of well tested 
tools for the computation of top, bottom and charm-pair hadropro duct ion at high 
energies, which has been developed for and used extensively in the LHC environment. 
For estimating lepton fluxes from atmospheric charm, the core of our calculation is 
an updated estimate of the NN ^ cc production cross-section in perturbative QCD, 
where N indicates a nucleon We use the QCD predictions for the inclusive cross- 
section pp —>■ cc at NNLO, and their comparison to those at NLO, to study the 
stability of the perturbative expansion as a guide for hxing parameters and inputs 
to be adopted in the application to atmospheric charm. In detail, this includes the 
PDF dependence, the choice of the central values for renormalization and factoriza¬ 
tion scales fiR and fip and the charm mass, as well as plausible intervals for their 
variations, given the fact that the global uncertainty bands at NLO due to scale, PDF 
and charm mass variation are large. Our hndings are summarized in the sequel. 

Fig. 3 displays the dependence of the total cross-section app^cc on the laboratory 
energy Eiab- The computation is performed in the hxed flavor number scheme (FFNS) 
with the number of flavors uj = 3, implying that charm is considered as a heavy 
state consistently included in the matrix elements with its mass different from zero 
and its presence excluded from the PDFs. The computation is performed in the 
theoretical framework as implemented in the HATHOR code [39]. Fig. 3 applies two 
different schemes for the heavy quark mass renormalization, the commonly chosen 
on-shell scheme with the pole mass and the MS scheme with the running mass 
rricifiR), where the renormalization scale for the evaluation of the mass has been hxed 

The experimental data in Fig. 3 for the hxed target experiments with energies 
up to Eiab = 10^ GeV are taken from Ref. [48] and for HERA-B from Ref. [49] (pur¬ 
ple points in Fig. 3). RHIC data from PHENIX and STAR have been published in 
Refs. [50, 51] (black points in Fig. 3) and LHC data are available from ALICE [11], 
ATLAS [12] and LHCb [13] (blue points in Fig. 3), see also Ref. [15]. Fig. 3 demon- 

^We use the approximation p ~ n ~ iV, neglecting mass differences between protons and neutrons 
by approximating all masses as nip. At high energies differences in the partonic content of p and 
n may also be safely neglected, whereas at low energies differences in the PDFs of these partons 
imply differences between pp and NN cross-sections up to a factor, depending on the specific PDF 
set, of a few percent at Eiab = 50 GeV, reducing to a few per mil above Eiab = 500 GeV. 
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Figure 3. Theoretical predictions for the total pp —?■ cc cross-section as a function of 
the laboratory energy Eiab at LO (dotted), NLO (dashed), NNLO (solid) QCD accuracy 
in the pole mass (left) and in the MS mass scheme (right) using the central set of the 
ABMll PDFs in the FFNS with Uf = 3. The scales were chosen as Pr = pr = 
with = 1.4 GeV in the on-shell scheme and as pR = pp = 2mc{mc) with 

mc{mc) = 1.27 GeV in the MS mass scheme, respectively. See text for details and 

references on the experimental data from fixed target experiments and colliders (STAR, 
PHENIX, ALICE, ATLAS, LHCb). 




Figure 4. Sensitivity of the total cross-section for pp — >■ cc to the factorization scale 
Pf at LO (dotted), NLO (dashed), NNLO (solid) QCD accuracy, in the pole mass (left) 
and in the MS mass scheme (right). The charm mass and PDFs were fixed as in Fig. 3. 
The central line at each order denotes the choice pR = pp- The upper and the lower lines 
at NNLO denote the cross-sections from the mass variation = 1.40 ± 0.15 GeV and 
mc{mc) = 1.27 ± 0.03 GeV, respectively. The arrows indicate the scale pR = pp equal to 
(left) and 2mc{mc) (right), respectively. 


strates the stability of the perturbative expansion of the Upp^cc cross-section through 
NNLO up to very high energies and good consistency of the predictions with the ex¬ 
perimental data. 
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Figure 5. Sensitivity of the total cross-section for pp —>■ cc to the factorization scale pp 
with the same PDFs and charm mass central values as in Fig. 4. The central line at each 
order denotes the choice pp = pp, the upper and the lower line the choices pp = pp/2 
and Pp = 2pp, respectively. The vertical bars give the size of the independent variation of 
Pp and pp in the standard range < pp, pp < 2mc°^'^ and mc{mc)/2 < pp, pp < 

2mc(m,c), respectively, with the restriction that 1/2 < pp/pp < 2. Again, the arrows 
indicate the scale pp = pp equal to 2mc°^'^ (left) and 2mc{mc) (right). 




Figure 6. Dependence of the total cross-section for pp —)• cc on the PDF choice at LO 
(dotted), NLO (dashed), NNLO (solid) QCD accuracy in the MS mass scheme. The charm 
mass and the scales were fixed as in Fig. 3. The upper and the lower lines at NLO and 
NNLO indicate the total Icr PDF uncertainty band for ABMll (left) and NNPDF3.0 PDF 
set (right) with rej = 3. Experimental data are the same as in Fig. 3. 


Related to the heavy quark mass renormalization is the choice of the numerical 
value for the charm quark mass. The Particle Data Group (PDG) [52] reports a very 
precise value of mc{mc) = 1.275 ± 0.025 GeV in the MS scheme. In case of charm, 
the conversion of the MS to the pole mass suffers from well-known convergence 
problems, see, e.g., Ref. [53]. In addition, the dehnition of the pole mass is based on 
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the unphysical idea of quarks as asymptotic states of the S'-matrix, so the accuracy 
in the pole mass is limited to be of the order (9 (Aqcd) by the renormalon ambiguity. 
The comparison of the Cpp^cc cross-sections in the two mass renormalization schemes 
at the nominal scales Hr = equal to and 2mc{mc) and taking the result 

with the precise PDG value as a reference, motivates our choice for the charm pole 
mass = 1.40 ± 0.15 GeV, as illustrated in Fig. 4. 

The behavior of the total cross-section as a function of the factorization scale 
/ii? is further explored in Fig. 5 by considering three different renormalization scales, 
Hr = Hf/‘2, Hf and 2hf-i at LO, NLO and NNLO QGD. Fig. 5 demonstrates, that the 
choice of the mass renormalization scheme is important, because the MS scheme leads 
to predictions with slightly improved convergence. Scale stability of the perturbative 
expansion at NNLO is reached in both schemes for scales /tr ~ /xj? > 2 GeV. As 
shown in Fig. 5 the use of the running mass scheme leads to a somewhat reduced 
scale uncertainty band at NNLO for the independent variation of hr nnd hf fh 
the standard range HR/Fncijnc) and HF/Fricijnc) G [1/2,2] and restricting the ratio 
1/2 < Hr/ f^F < 2, as compared to the pole mass scheme. Similar features, although 
much more pronounced, have been found already for the ti hadroproduction cross- 
sections and differential distributions in Ref. [40] . Both in the pole and in the running 
mass scheme the point of minimal sensitivity, i.e., the region where the cross-sections 
predictions at NLO and NNLO approximately coincide, turned out to be around 
scales Hr = f^F ^ and larger. This justihes scale choice adopted in Fig. 3. 
Translating this value into a dynamical scale, more suitable to describe the dynamics 
of heavy-quarks in differential distributions, we will use in the following a dynamical 
central scale for our calculation hxed to hf = HR — \/Pt + where pt is the 
transverse momentum of the emitted charm quark. 

Another important factor contributing to the theoretical uncertainties of the 
cc hadroproduction cross-section originates from the choice of the PDF set. We 
have taken predictions at NNLO accuracy as the basis of our central PDF choice. 
Among the different possibilities, currently available in the LHAPDF interface [54], 
we have chosen the ABMll one [55], together with the respective value for the 
strong coupling constant as{Mz), as a default. In the FFNS with Uf = 3, this 
PDF features a central set complemented by 28 variations, allowing to estimate a 
PDF uncertainty band at the la level. The ABMll PDFs are compatible with 
ABM12 [56], where the latter set has been tuned to LHG data. Moreover, the 
predictions of the ABMll and ABM12 family for gluon PDF at low Bjorken-x values 
are in complete compatibility with the only PDF £t available in literature so far 
including LHGb data on cc and bb hadroproduction, that has recently been performed 
by the PROSA collaboration [57, 58] 

As an alternative, we have used the 3-flavor central PDF set at NLO available 
®See also footnote 6. 
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from CTIO [59], which also provides results characterized by partial compatibili¬ 
ty with the PROSA £t (differences lie within about 2a) and with ABMll. At 
NNLO, both the CTIO and the ABMll PDF sets give positive results for the to¬ 
tal cross-section for cc hadroproduction in the highest energy range ranging up to 
Eiab ~ 10^° GeV, together with an uncertainty band for the PDFs which always stays 
positive as well, see Fig. 6 (left) for predictions from ABMll. 

On the other hand, PDF sets with unconstrained gluons at small x lead to 
very different results. In particular, the NNPDF3.0 set [60], characterized by a 
different parameterization, leads in case of the highest energies, to a huge uncertainty 
band, even covering a range with negative cross-section values. As shown in Fig. 6 
(right), the cross-sections obtained with NNPDF3.0 at NNLO do not remain positive 
anymore already for Eiab ^ 5 ■ 10^ GeV, an energy well below the one so far covered 
in run 1 by the LHG Thus, in the remainder of this paper we only consider the 
central value of the NNPDF3.0 set, with the purpose of quantifying differences with 
respect to the central values of the other families. 

The PDFs from MMHT [62], and in particular their central £t, lead to negative 
cc hadroproduction cross-sections at NNLO for energies above Eiab ~ 5 x 10® GeV. 
Such an unphysical feature also affects predictions for the longitudinal structure 
function El in deep inelastic scattering as noted, e.g., in Ref. [63]. While the MMHT 
set seems to be valid for the description of the production of heavier particles at 
LHG energies, its extrapolation to higher energies characterizing several astroparticle 
physics problems is quite questionable. Thus, we neglect this PDF family in the 
present study. 

Finally, the small x region is not only of importance for the gluon PDF, but 
also when considering the behavior of the perturbative hard parton scattering cross- 
section. The high-energy factorization of the cross-section [64, 65] in the limit when 
the center-of-mass energy S is much larger than the heavy-quark mass provides an 
effective theory for the description of the high-energy logarithms in S/m?. These 
behave as \v?{S/m?) ~ const, at NLO, as ln^(S'/m^) at NNLO, and so on, see, e.g.. 
Ref. [66] and studies of operator matrix elements in deep-inelastic scattering at three 
loop order in the small-x limit [67, 68]. At the energies currently considered, even up 
to Eiab 10^° GeV, their numerical importance is, however, strongly suppressed in 
the convolution integral of the hard partonic cross section at small x with the large 
X part of the gluon PDF (and vice versa). The apparent convergence of perturbative 
expansion for the app^cc cross-section through NNLO observed in Fig. 3 underpins 
this fact. 

® Data for the hadroproduction of heavy quarks at the LHC can therefore be used to further 
constrain these PDFs at small x. Very recently, following a research line non too different from 
the one already pointed out by the PROSA collaboration, Ref. [61] has considered constraints on 
the small-a; gluon from charm hadroproduction at the LHC, working in a scheme with 5 flavors, 
though. In contrast. Refs. [57, 58] and our work make consistent use of the FFNS with 3 flavors. 
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The cc differential distributions which we use in this paper are at parton level 
exact to NLO in QCD, because differential predictions for cc hadropro duct ion are not 
yet available at NNLO. For generating these distributions we use the POWHEG-BOX [46], 
complemented by the event generator PYTHIA-6.4.28 [14], in a pr-ordered tune 
belonging to the family of Perugia tunes [69], for describing parton shower and 
hadronization. This provides us with differential distributions of H-hadrons at 
NLO accuracy in QCD with NLO matching to parton showers according to the 
POWHEG method. Beyond the resummation to leading logarithmic accuracy provided 
by the parton showers approaches, next-to-leading logarithmic (NLL) contributions 
of as obtainable by an approach like FONLL, are not considered here. This 

is justified because the computation of Z-moments requires an integration over the 
whole kinematically accessible range in px and thus exhibits less sensitivity to the 
shape of the pt distribution at large pt, which is mainly influenced by the NLL 
corrections provided by the FONLL approach. 

Moreover, following Ref. [9], we derive the total charm cross-section in the ine¬ 
lastic proton-Air {pA) collisions, from the pN cross-section, by using the formula 


(2.9) 



where A = 14.5 for a nucleus of air. Here, we take 7 = 1 assuming a linear su¬ 
perposition since, for light nuclei, the effect of nuclear shadowing is expected to be 
small [15, 70]. 


2.4 Z-moments: Zph, Z^i, Zpp, Zhh 
2.4.1 Zph 


The inclusive cc cross-section is an essential ingredient for the estimate of the charm 
contribution to lepton fluxes. The latter ones, in fact, depend on the Z produc¬ 
tion moments which can be expressed as integrals over the differential distribution 
dap^charm{,E/ xe)/ dxEi through the formula. 



( 2 . 10 ) 


with the ratio xe = E/Ek. Here, E^ is the nucleon energy in the laboratory frame 
and E the energy of the produced particle (charm). The primary CR nucleon flux 
is (j)p and we have assumed that charm is all produced in cc pairs from primary CR 
nucleons (denoted by p) interacting with the Earth atmosphere, i.e.. 


dCp^charm/dXE 2 d(JpA-^cc /dXE , 


( 2 . 11 ) 


and that <Jp{E) in eq. (2.10) coincides with the total inelastic proton-Air cross- 
section crp7f4jr(E). The lower integration limit would ideally correspond to the case 





of Ek ^ + cx). We thus replace it with e, as we compute the cross-sections for 

Ek limited to a hnite value, with e small enough that the results for Z-moments are 
almost independent of its variations e^ar < e. 

As discussed in Sec. 2.3, for the differential cross-sections dapA-^cc/dxE and 
dapp^cc/dxE no predictions at NNLO are available at present. Thus we compute 
it at NLO through POWHEG-BOX, using PDFs and parameters such as rric and scale 
choices as described above, so that we are close to the point of minimal sensitivity 
for the total cross-sections, where the differences between NLO and NNLO QCD 
predictions are small. 

The hadronic moments Zph with h = Df and were calculated in 

Ref. [9] from the partonic moment Zp^c/iarm, through the relation Zp/j = fcharm,h Zp^charm, 
including conversion factors fcharm,h which represent the fraction of charm quarks 
emerging as specihc hadrons after fragmentation. In Ref. [15] a more rehned ap¬ 
proach was used: differential distributions for hadrons were obtained from those for 
quarks after convoluting the latter ones with fragmentation functions that were as¬ 
sumed to depend on the energies through the ratio E^/Echarm and to be independent 
of the beam energy. On the other hand, the use of POWHEG-BOX allows us to follow 
a different path, i.e., we take into account parton shower and fragmentation effects 
by means of the Monte Carlo PYTHIA event generator applied to the Les Houches 
events at hrst radiation emission level obtained by running POWHEG-BOX. This al¬ 
lows us to directly extract differential distributions dapp^h+x/dxE for D-hadrons 
{xe = Eh/Ek), whose shape may, in general, differ from those of the charm quarks 
dapp^c+x/dxE {xe = Echarm/Ek), implying that a global rescaling factor is too naive 
an approximation for the translation of quark distributions at parton level into the 
corresponding ones at hadron level. 

As an example, the differential distribution dapp^h+x/dxE is shown for the case 
of hadrons in Fig. 7 for a pp —)■ cc collision characterized by Eiab = 10^ GeV. 

The uncertainties due to scale and mass variation around the central predictions are 
shown in the left and right panels of the hgure, respectively. They are, in general, 
large. The scale uncertainties turn out to be almost constant with xe, whereas 
the uncertainties due to the variation of the mass increase with increasing xe on 
kinematical grounds. Differential distributions for different Eiab show a qualitatively 
similar shape, a scaling property already pointed out in the literature. 

Uncertainties due to the PDF variation are shown in the left panel of Fig. 8. For 
illustration, the predictions for the central ht and the 28 additional variations which 
parametrize the PDF and uncertainties of the ABMll PDFs with = 3 in the 
FFNS, are displayed individually 

^All differential results for PDF (scale) variations in this paper have been obtained after showe¬ 
ring with PYTHIA sets of Les Houches events generated in the POWHEG-BOX framework by explicitly 
fixing different PDFs (scales) in the input cards, without making use of reweighting facilities. This 
ensures a fully consistent computation of the Sudakov form-factors, including the specific PDF set 
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Figure 7. Differential distribution da/dxE for pp ^ cc ^ D^-\-X from POWHEG-BOX inter¬ 
faced to PYTHIA at Eiat = 10^ GeV. Central scales were fixed as pn = pp = ■\JPr^c + > 

central mass as = 1.4 GeV, and PDFs as the central set of the ABMll NLO family 
with Uf = 3. The uncertainty bands related to scale variation (at fixed and PDFs) 

and mass variation (at fixed pp, pp and PDFs) are shown in the left and right panels, 
respectively. The lower subpanels display the band for the relative uncertainties when 
normalized with respect to the central prediction. 

As is evident from Fig. 8, the differences with the central values of other PDF 
sets (CTIO and NNPDF3.0 at NLO in the rif = 3 FFNS with their respective default 
value for as{Mz)) turn out to be larger than the combined PDF and as uncertainty 
coming from the 28+1 ABM individual sets. This reflects a feature already observed 
at NLO for several examples of LHC cross-sections, where differences arising from the 
use of different PDF families turned out to be larger than those from the variation of 
as and the PDFs through the sets belonging to a same family, see, e.g.. Refs. [55, 71]. 

These differences propagate to the computation of the Z-moments, although the 
latter quantities are integrals over all possible xp values. As an example, in the right 
panel of Fig. 8, the Z-moments for hadroproduction are shown as a function 
of the energy in the laboratory frame for the different PDF choices just 

discussed above (left panel. Fig. 8). Here, the power-law spectrum as been chosen 
as input for the CR flux. Thus, the change of shape visible around 5 TO® GeV is 
due to the change of the spectral index in the power-law spectrum around the knee. 
The largest differences between different PDF sets appear at the lowest and at the 
highest D® energies. 

It is worth noting that pp collisions with Eiab > Eiab^po contribute to the Z- 
and scale in the whole integrand, in the separate generation of each set of events. 
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Figure 8. Differential distribution da/dxE for pp ^ cc ^ + X from POWHEG-BOX 

interfaced to PYTHIA at Eiab = 10^ GeV (left) and Z-moment for hadroproduction 
as a function of (right). Scales were fixed as pR = pp = mass as 

^poie — I ^ GeV. The red lines correspond to the central fit and the 28 additional sets 
for parametrization of the PDF and ag uncertainties of the ABMll PDFs at NLO with 
Hf = 3. Predictions using central fits of the CTIO and NNPDF3.0 sets at NLO, with 
Uf = 3 each, are also shown (in blue and light-blue, respectively), together with their ratio 
to the predictions with the central set of ABMll at NLO. 

moment at any given fixed energy Eiab^D^- Although, in line of principle, Eiab can be 
very large, in practice it turns out that the largest contribution comes from values of 
Eiab within the range Eiab^po < Eiab < (100 — 1000 ) x Eiab, do, due to the fact that the 
distribution in xp = Epo/Eiab is rapidly suppressed for large xp. As a consequence, 
for energies as those probed by IceCube, the contributions to the Z-moments come 
mainly from regions with a center-of-mass energy a/A non too high with respect to 
the energy range reached and probed so far at the LHC, where perturbative QCD 
in the standard formalism of collinear factorization has been tested to work. Any 
deviations from this formalism which may exist at the highest energies, e.g., in the 
form of non-linear effects (like gluon recombination as opposed to gluon splittings) 
or due to the dominance of large logarithms ln(S'/m^) subject to resummation on the 
basis of a different factorization formalism {kp factorization) [65], are, thus, expected 
to have only a small impact on the Z-moments we are interested in for the aim of 
understanding the IceCube results. 

The relative importance of Z-moments of different H-hadron species is shown 
in Fig. 9, where the Z-moments of positively charged H-hadrons and are shown. 
The contribution is the dominant one at all hadron energies. The different shape 
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Figure 9. Z-moments for the hadroproduction of selected L)-hadrons containing a 
c valence quark (electric charge Q = +2/3) {D^, D~^, and A+). Charm mass, {fiR, 
scales and PDFs were chosen as the central values in Figs. 7 and 8. The power-law CR 
spectrum has been used as input of our calculation. 


of the A+ contribution with respect to those of other hadrons is partly related to the 
shape of the differential distribution da/dxE of this hadron, which turns out to be 
less steep at large xe than those of the D-meson distributions. 

The expressions for Zph enter directly into those for the fluxes eqs. (2.4) and 
(2.5), obtained after solving the system of coupled differential equations describing 
the linear development of the hadronic cascade in the atmosphere, under the ap¬ 
proximations outlined in Sec. 2. 

2.4.2 Zhi, Zpp and Zhh 

In the following we briefly summarize our treatment of the other Z-moments Zhu 
Zpp and Zhh entering eqs. (2.4) and/or (2.5). 

For Zhi, our treatment of the semileptonic decay of D-hadrons follows closely 
Ref. [15]. Form factors for analytical decay distributions h —?■ were extracted 

from Ref. [7] and for the decay branching ratios the most recent values reported by 
the PDG [52] were taken. 
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In order to evaluate the proton regeneration Z-moment, Zpp, we have approxi¬ 
mated the inelastic xe distribution for the process pA —)■ pX by a scaling form 
da/dxE ~ — xe)"'{^ + n) with n = 0.51, as already done in Ref. [15]. 

Here, for we have considered the three different models already described in 
Sec. 2.2, which also enter the generation of the production moments Zph- 

Due to the obvious difficulties in measuring hA cross-sections with h being a D 
or B hadron, caused by the short lifetime of these particles, the moments Z^h are 
approximated by considering the available estimates for KA cross-sections, on the 
basis of analogies between K and D mesons. Both include quarks belonging to the 
same flavor family (charm quarks in case of D’s are replaced by strange quarks in case 
of i^’s). In particular, following Ref. [15], the inelastic xe distribution for the reaction 
hA —)■ hX is estimated by the ansatz da/dxE ~ a)!^^{Eiab) (1 — xe)"' (1 + n), 

where n = 1 and cr^^N total inelastic cross-section for R'^-nucleon interactions. 
To estimate the latter one, total and elastic cross-sections were extracted from the 
latest version of the PDG. However, the behavior of the K^p elastic cross-section at 
high energies is uncertain because no data above Eiab ~ 300 GeV exist. We have thus 
assumed that the slope of the K^p elastic cross-section at high energies is similar 
to the one of the pp elastic cross-section, which was recently constrained at LHG 
energies by TOTEM data [72]. 

The regeneration Z-moments Zpp and Z^^ enter the expression for the atte¬ 
nuation lengths Ap and A/^, respectively, as dehned below eq. (2.5). It is worth noting 
that estimates of the Z/i/j-moments and of the related uncertainties only affect the 
high-energy approximated solution of the cascade equations, eq. (2.5). 

3 Neutrino fluxes and their uncertainties 

The IceGube experiment is looking for a diffuse flux of neutrinos at high-energies, 
including both downward and upward going neutrinos, by trying to establish the 
nature of the observed events as due to astrophysical signals or to the atmospheric 
(conventional -|- prompt) background. In this Section we focus on (i/^ -f- i/p)-fluxes, 
by taking into account that the largest contribution to the atmospheric conventional 
neutrino flux (to which we will compare our prompt flux) comes from this flavor. 
Predictions for other leptons can be obtained with the same method as well. The 
qualitative/quantitative difference between the results for + Vp) fluxes and those 
for other leptons depends on the specihc decay modes and branching fractions of D 
hadrons in each species. 

The lepton fluxes are derived after evaluating all quantities entering eqs. (2.4) 
and (2.5), already described in previous Sections, and by interpolating between the 
high energy and the low energy solutions according to eq. (2.6). An example of the 
typical behavior of the two solutions and of their interpolation is shown in Fig. 10 
for the case of a power-law primary GR spectrum as input of the whole calculation. 
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vl^ + anti-v^ flux 



Figure 10. The + P^)-flux as a function of the neutrino energy Ei^b^u illustrating 
the geometric interpolation between low-energy and high-energy solution to the cascade 
equations, in case of a power-law primary cosmic-ray spectrum. 

In the following we will present the central values of our fluxes, together with the 
uncertainty bands arising from the different source of uncertainties, both of QCD 
and of astrophysical origin. 

3.1 Main uncertainties from QCD and astrophysics 

Uncertainties on the fluxes whose origin can be ascribed to perturbative QCD mainly 
reflect those uncertainties already found in the differential distributions da/dxE and 
in the Z-moments. In particular, we discuss in the following the scale, charm mass 
and PDF variation, as well as matching uncertainties, related to the NLO matching 
to the parton shower. 

Uncertainties in the + P^)-fluxes due to and scale variation for a 
power-law CR spectrum as input, are reported in Fig. 11, while the corresponding 
uncertainties for other CR spectra are shown in Fig. 12. The scale variation turns 
out to be the largest source of uncertainties. Including cases with fip ^ fip, i-e., the 
independent variation of fip and fip, leads to an uncertainty band which is almost 
uniform on a wide interval of energies Eiab- In this respect our hndings in Figs. 11 
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Figure 11. The + ^'^)-fluxes as a function of the neutrino energy Eiab^u with uncer¬ 
tainties due to renormalization and factorization scale variation. Charm mass, PDF and 
scales were fixed as in the left panel of Fig. 7. 

and 12 are different from the result of Ref. [15], where the non-diagonal choices with 
hR 7^ hF were neglected and the scale uncertainty is underestimated, especially at 
low energies. 

Uncertainties arising from the variation of the charm mass rric within the range 
motivated in Sec. 2.3 are illustrated in Fig. 13 for a power-law CR spectrum and 
in Fig. 14 for other CR spectra. The mass variation turns out to be the second 
largest source of QCD uncertainties, with an uncertainty band slightly decreasing 
for increasing laboratory energies Eiab- 

Uncertainties in the neutrino fluxes related to the PDF variation are displayed 
in Figs. 15 and 16 in case of a power-law primary CR spectrum and for the different 
variants of Gaisser spectra, respectively. As already discussed for the case of the 
Zp/j-moments, the difference of the predictions with the ABMll set (3-flavor FFNS) 
and the central value of other PDF families (CTIO and NNPDF3.0 at NLO with 
Uf = 3) turn out to be larger than those coming from the 28 sets in the ABMll 
£t for the combined PDF and as uncertainty. While the neutrino fluxes from the 
different PDF fits look quite consistent among each other for energies in the interval 
10^ < Eiab, u < 4:- 10“^ GeV, visible differences between the fluxes from different PDF 
families start to appear at higher energies. In this region, the predictions based on 
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Figure 12. Same as in Fig. 11 for different primary CR spectra, where each panel 
corresponds to a variant of the Gaisser primary spectrum, cf. Sec. 2.1. 

the ABMll PDFs are the smallest ones, which is related to differences in the shape 
of the gluon PDF, the nominal values of the strong coupling as{Mz) at NLO being 
largely the same among the ABMll, CTIO and NNPDF3.0 sets used in this work. 

Finally, we provide a hrst estimate of the uncertainties in the NLO matching to 
the parton shower (NLO + PS) by varying the hdamp value in the POWHEG-BOX, which 
parameterizes the freedom in choosing the form of the separation of the NLO real con¬ 
tribution R into a singular piece plus a piece damped in the singular region and thus 
treatable as a hnite remainder [73], R = Rg + Rf, with Rg = R h'^amp / i^damp + Pt) 
and Rf = R / {h?damp + Pt) Only Rg enters the exponent of the Sudakov 
form factor and hdamp = +oo corresponds to the default choice in POWHEG-BOX so that 
R = Rg, whereas the limit hdamp —t 0 allows to decrease the amount of radiation that 
is exponentiated and to recover the dependence (pure NLO) in the high-pj’ limit. 
In this work we use variations of hdamp in the interval {rric, 2mc, 4:mc, +C)o}. This is 
inspired by similar choices performed in experimental studies of ti hadroproduction, 
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Figure 13. The + ^'^)-fluxes as a function of the neutrino energy Eiab^u with uncer¬ 
tainties due to the variation of the pole mass = 1.40 ± 0.15. PDF, scales, 

charm mass were fixed as in the right panel of Fig. 7. 

see, e.g., the ATLAS note [75]. The uncertainty is estimated as the envelope of the 
predictions corresponding to the different choices above, as shown in Fig. 17 in case 
of a power-law primary CR spectrum and in Fig. 18 in case of the Gaisser spectra. 
Although several discussions are still on-going about the most meaningful way of 
providing NLO -|- PS matching uncertainties which is why we consider our result 
as a first rough estimate, we would like to point out that the uncertainty we got is 
quite small (less than 10%) with respect to other uncertainties of QCD origin. This 
is related to the fact that the key quantities in perturbative QCD to compute Z- 
moments, the differential cross sections da/dxE, are integrated over the entire range 
of transverse momenta pt- We thus believe that this conclusion is quite robust, i.e., 
it does not depend on very specific details of the way the matching uncertainty is 
estimated. 

A summary of the main QCD uncertainties (mass, scale, PDF variation) in re¬ 
lation to uncertainties of astrophysical origin, in particular those arising from the 
variations of the primary CR flux used as input in the + P^)-flux calculation, is 
provided in Fig. 19. In the first three panels of this figure we show separately the 
uncertainties due to scale, mass and PDF variation (by restricting ourselves to the 
ABMll PDF set) for all hve primary CR fluxes considered as input in this paper. 


23 















+ anti-v^ flux 


Vp + anti-Vp flux 




v^ + anti-Vp flux 


+ anti-v^ flux 



Eiab,v. (GeV) 


i^charm - f -40 GoV 
i^charm ~ f -25 GeV 
*^eharm ~ f GeV 



10^ 10^ 10® 10® 10^ 10® 
E|ab,u ( GeV) 



Figure 14. Same as in Fig. 13 for different primary CR spectra, where each panel 
corresponds to a variant of the Gaisser primary spectrum, cf. Sec. 2.1. 

The uncertainties due to scale variation are the dominant component. Apart from 
that, it is evident that at energies > 10® GeV uncertainties due to variations in 
the CR fluxes dominate over those from mass and PDFs. On the other hand, un¬ 
certainties related to QCD effects always dominate at energies < 10® GeV, where 
the primary CR fluxes are well constrained by several measurements (see Sec. 2.1). 
Finally, in the last panel of Fig. 19 we show the quadrati c combination of the un- 
certainties above, assumed as independent, i.e., Aqcd = \J^‘pdf- 
For Eiab,u = 10® GeV, —72% < Aqqd < +84%, i.e., the uncertainty is slightly 
asymmetric, and it slightly changes (a few percent) at higher energies. 

We also observe that with a restricted scale variation interval, neglecting the com¬ 
binations {fiR, fir) = (0.5, 1) and (1, 0.5) po for po = \Jpt,c + 4m2, scale uncertainties 
and, as a consequence, the total ones, are reduced, as shown in Fig. 20. In particular, 
for Eiab,u = 10® GeV, the combined uncertainty amounts to —48% < Aqcd A +63%, 
and it changes by a few percent at higher energies. 
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Figure 15. The + i/^)-fluxes as a function of the neutrino energy Eiab,v with un¬ 
certainties due to PDF variation in the 3-flavor ABMll PDF set at NLO (red band) 
and predictions for the central set of the 3-flavor NNPDF3.0 (light-blue line) and CTIO 
(solid blue line) PDFs at NLO, respectively. Charm mass and scales were fixed 

as in Fig. 8. The power-law cosmic ray flux has been used as input in the calculation of 
Z-moments. 

3.2 Other uncertainties 

In the previous Sec. 3.1 we have provided a minimal estimate of the combined QCD 
and astrophysical uncertainties which affect our predictions for -f- P^)-flux. In 
the following we shortly describe other sources of uncertainties which could be added 
to the previous ones. 

A further QCD contribution arises from heavier hadrons, in particular S-hadrons, 
which are also a source of neutrinos, but whose effect has been neglected here. Given, 
that the cross-section for bb hadroproduction with respect to the one for cc hadropro- 
duction is smaller by a factor of order 20 at LHC energies and still suppressed by 
a factor of order 10 at Eiat = 100 TeV, we expect that the bottom-quark contribu¬ 
tion can be neglected with respect to the charm one at the energies of interest for 
IceCube. However, bb hadroproduction may play a larger role at ultra-high-energies. 

Other uncertainties can be attributed to the approximate description of the de¬ 
cay of heavy hadrons. In particular, a component of secondary neutrinos coming 
from the decay of the lighter mesons (baryons) produced as decay products of D- 
mesons (baryons), is missing in our computation as well as in many previous ones 
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Figure 16. Same as in Fig. 15 for different primary CR spectra, where each panel 
corresponds to a variant of the Gaisser primary spectrum, cf. Sec. 2.1. 

(see e.g. Ref. [15]). Furthermore, from the QCD point of view, non-perturbative ef¬ 
fects, suppressed by powers of Aqcd/i^, are increasingly important for smaller quark 
masses m. In this respect, one should account for a contribution to the uncertainties 
due to fragmentation, arising from both fragmentation fractions and fragmentation 
functions [76]. The latter can be estimated, for instance, by varying the choice of 
the functional form of fragmentation functions for heavy flavors in PYTHIA together 
with the parameters involved. The corresponding uncertainty estimates have been 
discussed in the literature, see, e.g., Ref. [41]. Potential uncertainties related to 
the variation of the partonic intrinsic transverse momentum (/ct) ~ Aqcd are in¬ 
deed smaller since they mostly affect the pt distributions to which our work is not 
particularly sensitive. 

From an astrophysical point-of-view, further uncertainties to be included encom¬ 
pass those related to a change in the p-Air cross-section, which affect both the Zph 
production moments and the Zpp regeneration moments, and, as a consequence, the 
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Figure 17. The + P^)-fluxes as a function of the neutrino energy Eiab^u with the un¬ 
certainties from the NLO -|- PS matching estimated in the POWHEG-BOX framework through 
hdamp variation. See text for more detail. The power-law cosmic ray flux has been used as 
input in the calculation of Z-moments. 

lepton fluxes. To that end, we consider the theoretical predictions coming from the 
three different models described in Sec. 2.2 (QGS JetO. Ic, SYBILL2.1 and the analyti¬ 
cal model eq. (2.8)). The larger the inelastic cross section cr*"'®*(p-Air) is, the smaller 
are the predictions for our fluxes, as is evident when comparing Fig. 21 with Fig. 2. 
However, these global effects on neutrino fluxes turn out to be not too relevant, i.e., 
the uncertainties coming from the use of different models, amount to less than 10% 
over the whole Eiab, v energy range considered. They are therefore much smaller than 
those from QCD effects discussed previously and those from the choice of different 
primary CR spectra, as shown in Fig. 22. 

4 Comparison with previous results and astrophysical im¬ 
plications 

We compare our prompt (z/^ -|- P^)-flux with previous results in the literature. In 
particular, it turns out that our central fluxes are in between central predictions 
recently obtained by another group using the standard hard-scattering formalism in 
QCD [15] and older predictions provided in Ref. [10] by making use of the dipole 
picture. In particular, our central values are a few ten percent larger (~ 40% at the 
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Figure 18. Same as in Fig. 17 for different primary CR spectra, where each panel 
corresponds to a variant of the Gaisser primary spectrum, cf. Sec. 2.1. 

energies of interest for the IceCnbe experiment) than the central valnes in Ref. [15], 
that lie in any case within our quoted uncertainty band, for the various primary CR 
spectra already considered in that work, as shown in Figs. 23 and 24. Note, that 
Ref. [15] is based on a completely independent QCD calculation and on different 
inputs and methods. 

On the other hand, differences with older calculations, like the one in Ref. [8] on 
the basis of PYTHIA and including QCD hard-scattering effects at leading order only, 
are obvious, especially regarding the shape of the distributions. This is the case not 
only for the lepton fluxes, but already for the Z-moments for D-hadron production. 

Interestingly, our results are very well compatible with those from the dipole 
model of Ref. [10], altough the latter were computed with older sets of PDFs: as 
shown in Fig. 23, for Eiab,u > 10^ GeV, our central predictions are included in the 
uncertainty band of the latter, whereas the central predictions from the dipole model 
are included in our uncertainty band. 
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Figure 19. Summary of the main QCD and astrophysical uncertainties affecting our 
central predictions for + r'^)-fluxes. Uncertainties due to scale, mass and PDF vari¬ 
ation (considering the ABMll PDF and a* uncertainty band), are shown separately and 
combined for each of the five primary CR spectra, cf. Sec. 2.1. 


In order to infer a value for the transition energy Etrans where the prompt neu¬ 
trino flux overcomes the conventional one, in Fig. 24 we compare our prompt lepton 
flux with the conventional neutrino flux originally computed in Ref. [77] for a power- 
law CR primary spectrum and rescaled to one variant of the Gaisser spectra in 
Ref. [15]. We obtain Etrans = ■10® GeV. Interestingly, the central value lies 

well within the interval (4 -10® — 10®) GeV where the IceGube experiment did not 
observe any event after the full 988-day analysis [3]. In fact, the IceGube collabora¬ 
tion has reported an excess of neutrinos in the diffuse flux, all lying in the neutrino 
energy regions [0.3 — 4] 10® GeV and [1 — 2] 10® GeV. According to our predictions, 
the “empty” region of IceGube corresponds to the “conventional-prompt” transition 
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Figure 20. Same as panel 1 and panel 4 of Fig. 19, but for a restricted choice of ///j and 
variations. Here the scale uncertainty is obtained as the envelope of the combinations 
{fiR, fip) = (0.5, 0.5), (1, 1), (2, 2), (1, 2) and (2,1) ^q, disregarding the cases with 
{hr, hf) = (0.5, 1) and (1, 0.5) ho- 


region, i.e., the region where the contributions of conventional neutrinos and prompt 
neutrinos to the total neutrino flux become of the same order of magnitude. We 
thus believe that the “empty” region seen by IceCube so far, should not be empty, 
but actually dominated by prompt neutrinos. However, the IceCube error bars in 
the “empty” region are still quite large, and we stress that the accumulation of more 
statistics is necessary before judgment can be made, whether this lack of signal is 
just an artifact due to poor statistics or due to some other technical issue, or instead 
has a real physical interpretation. 

At higher energies, on the other hand, the total observed neutrino flux 4>{E^) 
for El, in the [1 - 2] 10® GeV energy interval looks to be slightly suppressed with 
respect to that in the [2 - 3] 10® GeV bins. However, looking at our central prompt 
flux distributions and summing them with the distributions for the conventional flux, 
as a hrst rough estimate it turns out that we would expect a much larger suppression 
in the [1 - 2] 10® GeV region with respect to the [2 - 3] 10® GeV one, disfavoring 
the interpretation that the events seen by IceGube in the [1 - 2] 10® GeV window 
are just due to a prompt neutrino component The difference between the IceGube 
(signal + background) observed total yield at high-energy and the yield for prompt 
neutrinos as predicted by our calculation is slightly reduced if we observe that our 
predictions have a sizable uncertainty band, meaning that even the shapes of the 
distributions can change in a non-negligible way when a higher-order calculation in 

®This interpretation is also disfavoured by IceCube observations of the arrival directions of the 
events with E > 6-10^ GeV, in presence of a /i veto (see Fig. 3 of Ref. [3]). 
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Figure 21. The + i>^)-fluxes as a function of the neutrino energy Eiab,v obtained by 
considering different models for the p-Air total inelastic cross-section shown in Fig. 2 for a 
power-law primary CR flux. Charm mass, PDFs and scale were fixed to our central values 
(see Figs. 7 and 8). 


QCD will be available, if we consider neutrino flux values corresponding to the upper 
value of our uncertainty band and if we use as input primary CR fluxes including 
a population of extra-galactic protons with very-high rigidity (i.e. variants 2 of 
Gaisser spectra, instead of variants 1 which have a mixed extra-galactic component 
with a lower global rigidity for the extra-galactic population). In fact, the latter give 
rise to neutrino spectra which are less severely suppressed at the highest energies 
than those from models with extragalactic mixed components, as is evident when 
comparing, e.g., the left and right panels of Fig. 24, obtained with the variants 1 
and 2 of the Gaisser 2014 spectrum, respectively. In order to go beyond these purely 
qualitative considerations and to draw more dehnite quantitative conclusions, one 
should dehnitely wait for more experimental statistics and, also insert our fluxes into 
the specific experimental analysis software. 

In any case, we would like to emphasize that the transition region for the prompt 
+ r';i)-flux in our calculation turns out to be also a transition region for uncertain¬ 
ties, i.e., the QGD uncertainties dominate the total uncertainties at energies below the 
transition region whereas the astrophysical ones start to give a progressively sizable 
contribution above it, pointing out the importance and necessity of pursuing further 
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Figure 22. Same as in Fig. 21 for the five different primary CR spectra considered in 
this work, cf. Sec. 2.1. 


studies of cosmic ray composition at the highest energies [78], and, possibly, future 
measurements independent of Monte Carlo simulations of hadronic-interactions at 
the highest energies. 

5 Conclusions 

We have computed the prompt neutrino fluxes from atmospheric charm using up-to- 
date theoretical results and tools for charm hadroproduction in perturbative QCD. 
Our results for the neutrino fluxes are several tens percent larger over a wide range of 
neutrino energies than predictions in the recent literature making use of Z-moments 
computed with the standard QCD hard-scattering formalism, that we also adopt. 
At the energies of interest for the IceCube experiment the increase of our prompt 
+ zx^)-flux amounts to some 40%. However, our uncertainties on the fluxes both 
of QCD and astrophysical origin are dramatically larger. Partly as an effect of this 
fact, even predictions obtained by making use of the dipole picture, representing an 
alternative description to the undelying hard-scattering, he within our uncertainty 
band over a wide energy range. 
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Figure 23. Comparison between the prompt + r'^)-flux obtained in this work (blue 
solid line with blue uncertainty band) with the central values of those previously obtained 
by other authors, for a power-law primary cosmic ray spectrum. The TIG flux (Ref. [8]) 
is shown by open magenta squares, the ERS central flux (Ref. [10]) and its uncertainty 
is shown in yellow, whereas the more recent BERSS flux (Ref. [15]) is shown by filled 
light-blue squares. 


We have discussed extensively the different sources of uncertainties which affect 
the fluxes. The main sources come from (i) the renormalization and factorization 
scale variation allowing for independent variations of ^ (ii) the charm mass 
uncertainties for the pole mass choice, and (iii) PDF uncertainties evaluated for the 
ABMll set and studied by comparing its predictions to the central predictions of 
different PDF sets (CTIO, ABMll, NNPDF3.0) at NLO. Further uncertainties due 
to hadronization and hadron decay have been discussed as well. In particular (i) and 
(ii) had not been included in a systematic way in studies in literature before, so we 
conclude that previous uncertainties on prompt neutrino fluxes are underestimated. 

The uncertainties of QCD origin dominate at low neutrino energies, whereas 
for increasing energies Eiab,u ~ 10® — 10® GeV the uncertainties in the astrophysical 
input, in particular the primary CR flux and its composition in terms of different 
populations, turn out to add a progressively important contribution to those from 
QCD. 

The results presented may beneht from a number of future developments. On 
the QCD side, a fully differential NNLO computation of charm hadroproduction. 
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Figure 24. Comparison between the prompt + r'^)-flux obtained in this work (blue 
solid line with blue uncertainty band) with the central values of the more recent BERSS flux 
(Ref. [15]) (light-blue squares) in case of recent primary cosmic ray spectra (Gaisser-2014- 
variant 1 on the left and Gaisser-2014-variant 2 on the right). The conventional neutrino 
flux computed by Honda (Ref. [77]), after its rescaling to the Gaisser-2014-variant 1 GR 
primary spectrum as presented in Ref [15], is shown by open circles. 

when available, will be of great help in reducing the theoretical uncertainties from 
scale, mass and PDF variation. In this respect, the role of resummation of different 
kinds of logarithms deserves further exploration as well. Furthermore, a dedicated 
systematic survey of the uncertainties related to both the fragmentation functions in 
the Monte Carlo parton shower matched to NLO predictions and the fragmentation 
fractions, would allow to quantify those effects. This could be a step towards the 
optimization of Monte Carlo tunes, to make them especially tailored to studies like 
those performed in this work. This optimization also concerns the search for the best 
parameter values for the description of dual and multiple particle interactions. On 
the experimental side, measurements of the cc and bb production cross-section at the 
LHC, looking not only at central rapidities but also in the forward rapidity regions, 
can be of importance especially at the highest energies, where the contribution of 
low X events becomes increasingly important. 

Finally, from the astrophysical point of view, one could obtain a substantial 
reduction of the uncertainties on prompt neutrino fluxes at the highest energies once 
issues related to the transition between a galactic and an extragalactic component 
in the CR primary spectrum and to the composition of the latter will be understood 
better. 

Our lepton fluxes will be made available as numerical tables for download at 
http: //www. desy. de/~promptf luxes. Further predictions can be requested to the 
authors of this paper by e-mail. 
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